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Abstract 

A new atsp2K module is presented for evaluating the electron density function of any multi- 
configuration Hartree-Fock or configuration interaction wave function in the non relativistic or 
relativistic Breit-Pauli approximation. It is first stressed that the density function is not a priori 
spherically symmetric in the general open shell case. Ways of building it as a spherical sym- 
metric function are discussed, from which the radial electron density function emerges. This 
function is written in second quantized coupled tensorial form for exploring the atomic spherical 
symmetry. The calculation of its expectation value is performed using the angular momentum 
theory in orbital, spin, and quasispin spaces, adopting a generalized graphical technique. The 
natural orbitals are evaluated from the diagonalization of the density matrix. 
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Operating systems or monitors under which the present version has been tested: HP XC System Software 
3.2.1, which is a Linux distribution compatible with Red Hat Enterprise Advanced Server. 

Programming language used in the present version: FORTRAN 90 

RAM: ?? MB or more 

Peripherals used: terminal, disk 

No. of bits in a word: 32 

No. of processors used: 1 

Has the code been vectorised or parallelized? : no 

No. of bytes in distributed program, including test data, etc.: ?? bytes 

Distribution format: gzipped compressed tar file 

CPC Program Library subprograms used: libraries of atsp2K 

Nature of physical problem 

This program determines the atomic electronic density in the MCHF (LS ) or Breit-Pauli (LS J) approxima- 
tion. It also evaluates the natural orbitals by diagonalizing the density matrix. 

Method of solution 

Building the density operator using second quantization - Spherical symmetry averaging - Evaluating the 
matrix elements of the one-body excitation operators in the configuration state function (CSF) space using 
the angular momentum theory in orbital, spin, and quasispin spaces. 

Restrictions on the complexity of the problem 

Original restrictions from atsp2K package, i.e. all orbitals within a wave function expansion are assumed 
to be orthonormal. Configuration states are restricted to at most eight subshells in addition to the closed 
shells common to all configuration states. The maximum size of the working arrays, related to the number 
of CSFs and active orbitals, is limited by the available memory and disk space. 

Typical running time 

The calculation of the electron density for a n = 9 complete active space (CAS) MCHF wave function 
(271 733 CSFs - 45 orbitals) takes around 9 minutes on one AMD Opteron dual-core @ 2.4 GHz CPU. 

Unusual features of the program 

The programming style is essentially F77 with extensions for the POINTER data type and associated mem- 
ory allocation. These have been available on workstations for more than a decade, but their implementa- 
tions are compiler dependent. The present code has been installed and tested extensively using the Portland 
Group, pgf90, compiler. 
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1. Introduction 



In electronic structure theory there are several approaches to describe the behavior of elec- 
trons in atoms and molecules. Most of them are based on the wave nature of the particles, 
permitting the system to be described by wave functions, as eigenstates of the Schrodinger equa- 
tion. The Hohenberg-Kohn (HK) 1 1|] theorems on the other hand say that the electronic structure 
of a system is completely determined by its ground state electron density function. According to 
the HK theorems, the energy of any system can be written as a functional of this density func- 
tion. Based on these results, within Density Function Theory (DFT), several methods have been 
developed to describe atoms and molecules through their density function [2]. The development 
of density functionals which yield a system's energy has become a major field of research in 
Chemistry and Physics. Nowadays a lot of research is being done to investigate how the elec- 
tron density function describes the system. In conceptual DFT for example, chemical reactivity 
indices are defined, which indicate how a system behaves in a chemical reaction, by considering 
perturbations to the electron density function. Although wave function methods were well de- 
veloped before, DFT is now the most widely used electronic structure method. The wide spread 
use of DFT can be accounted to the relative computational ease with which energies can be de- 
termined. Where a wave function describing an -particle system involves the position- and 
spin- coordinates of all electrons, a density function, describing the same system, only depends 
on the coordinates of one particle. Following the work of McWeeny |3J|, one can try to extract 
physically essential features from the electron density function. 

Some of the present authors have established the periodicity of the atoms in Mendeleev's 
periodic Table by making an information theoretical analysis of the electron density functions 
as probability distributions fl. Another work quantifies the relativistic effects on the basis of 
a comparison of density functions calculated within the one-configuration Hartree-Fock and 
Dirac-Fock approximations Ol • 

The present code is an extension of the Atomic Structure Package atsp2K | (]] for evaluating 
the atomic density function from non relativistic and relativistic (in the Breit-Pauli approxima- 
tion) multiconfiguration ab initio wavefunctions of atomic systems, adopting an efficient ap- 
proach for spin-angular integrations |3. 8]. It allows the investigation of correlation effects on 
the density function for any non-relativistic correlation model, and of relativistic effects in the 
Breit-Pauli approximation. 



In quantum chemistry, the natural orbitals (NO) are known to provide a particularly efficient 
choice of single-particle states JsKjil. Moreover, NO give the most rapidly convergent approx- 
imation to the total wave function and are often used as a basis set for generating a better wave 
function in an iterative manner. In atomic physics, NO are rarely used, although they consti- 
tute the orbital basis of the reduced form of the MCHF expansions for helium-like and nominal 
two-electron atomic systems fill . It would be worthwhile to study their potential for more than 
two electrons in the search of efficient optimization strategies. The present code fills this gap by 
building the natural orbitals through the diagonalization of the density matrix. 

2. On the symmetry of the density function 

In this section we start by formulating the multiconfiguration wave function for a well defined 
atomic state, and we calculate the corresponding density function. From this calculation, we 
regain the specific angular (non-spherical) dependence of the density function. We also present 
different ways for deriving a spherical electron density function. 

2.1. The multiconfiguration many-electron wavefunction 

In the multiconfiguration approach, the A^-electron wavefunction m l m s is a linear combi- 
nation of M configuration state functions (CSFs) ^> a ,LSM L M s which are eigenfunctions of the total 
angular momentum L 2 , the spin momentum 5 2 and their projections L z and S z , with eigenvalues 
h 2 L(L + 1) , h 2 S(S + 1), KM L and HM S , respectively 



The set of variables {x ; } represent the electron's space and spin coordinates Xj = (rj, crj) = 
(rj, §j, ipj, crj). The individual CSFs are built from a set of one-electron spin-orbitals, 



where R„i(r) = P„i(r)/r, Yi mi (ft,<p) and^ smi (cr) are the radial, the angular and the spin parts of 
the one electron functions. The mixing coefficients {c,} and the radial functions [R n i,(r)} are 
solutions of the multiconfiguration Hartree-Fock method in the MCHF approach. For a given set 
of orbitals, the mixing coefficient may also be the solution of the configuration interaction (CI) 



M 




(1) 




(2) 
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problem. The relativistic corrections can be taken into account by diagonalizing the Breit-Pauli 
Hamiltonian Jl^j in the LS /-coupled CSF basis to get the intermediate coupling eigenvectors 

M' 

"F tt /Af(xi. ■ • -XaO = J] Gi ®( a i L i s i JM > x u ■ ■ -x N ) . (3) 

i=l 

2.2. The non-spherical density function 

The so-called "generalized density function" |3 |or the "first order reduced density matrix" 
I13I is a special case of the reduced density matrix I I0II3I1 

yi(Xi,x'i) = N I ¥(Xi,X2, . ..,Xjyr)T*(x'i,X 2 , ...,x N )dx 2 ...dx N , (4) 



where T^xi , X2, ■ • ■ , Xjv) is the total wave function of an N electron system and T^Xi , X2, . . . , x#) 
is its complex conjugate. The spin-less total electron density function p(r) is defined as the first 
order reduced density matrix, integrated over the spin and evaluated for X] = x' 1 

/o(ri) = J~yi(xi,xi>fo-i. (5) 

This electron density function is normalized to the number of electrons of the system 

J p(r) dr = j 'p(r) r 2 sin &drd&d<p = N . (6) 

As discussed in IU3I1 . the single particle density function can be calculated by evaluating the 
expectation value of the 6(r) operator, 

p(r) = J~¥(xi,X2,...,x A r)(5(r) , F*(xi,X2, . . . ,x N ) dxrfx 2 . . .dx N , (7) 

where <5(r) probes the presence of electrons at a particular point in space and can be written as 
the one-electron first-quantization operator 



N 



6(f) = Yj 6(r ~ rd ■ (8) 



Expressing each S(r - r,) term in spherical coordinates [ 14] 



6(r - r,) = -5— - 6(r - n ) 6(§ - §d 6(<p - tpd , (9) 
r L sin 17 



and introducing the closure relation 



<fi)Y; n (&', ip') = 6(cos & - cos 6(<p - <p') , (10) 
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the operator ® becomes 



6(r) = £ <5(r - n) = Yj S( ~ r ~ r ') Z Ylm( ^ f )Y i'n (0 - <Pd ■ (ID 

The exact spin-less total electron density function (0 evaluated for an eigenstate with well- 
defined quantum numbers (LS MiM$), is 

1 N 

p{y) lsm l m s = ^ y imW _ <^ L5MiMs I Yj - n) rjtfu <Pi)\V aI sM L M s > (12) 

lm i= 1 

It is important to realize that the spherical harmonic components are limited to the /-even con- 
tributions, since the bra and ket states have the same parity n = (— 1)^' \ Applying the Wigner- 
Eckart theorem gives 



21 



p{r) LSM L M s = £ Ylo(§t(p) l_ 

/„„„=() r 



L I 
-M L M L 



(V*ls II J] S(r - n) Yj&u <pdW<*s ) 



(13) 



where 



/ \LSMlM* - i 1 1 



1 



L—Mr 



L 21 L 
-M L M L 



(ToM 1| £ <J(r - r,) Y* 2l (& u tpiWaU > ■ (14) 



This result recovers Fertig and Kohn's analysis [17] for the density corresponding to a well- 
defined (LS MiMs ) eigenstate of the Schrodinger equation. In this paper, the authors observed 
that the self-consistent field densities obtained via the Hartree and Hartree-Fock methods gener- 
ally violate the specific finite spherical harmonic content of p(r) LiSMiMs , They also mention that 
this exact form can be obtained by spherically averaging the effective potential, yielding single- 
particle states with good angular momentum quantum numbers. The atomic structure software 
pac-e^Kfl applied approach , as was done in the original atomic Hartree-Fock the- 
ory 11 (Si HQ . This implies two things: i) the density function p(r) iSMtMs calculated from 
any multiconfiguration wave function of the form ([1), is not a priori spherically symmetric, ii) 
this density function will contain all spherical harmonic components (up to 2L) as long as the 
one-electron orbital active set spanning the configuration space is /-rich enough. 



The same result can be obtained by reducing the many-electron reduced matrix element as a sum over one-electron 
reduced matrix elements as done in flql . 
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The density function can also be expressed in second quantization [2I]. Introducing the nota- 
tion q = n q l q m\m Sq for spin-orbitals, expression © becomes 



ri(xi,x'0 = Yj D pi ^(x'i)^ 9 (xi) , (15) 

pq 

where D pq are elements of the density matrix which are given by 

D pq = {^\a\a q \Y) . (16) 

The sum in eq. ([T5T > runs over all possible pairs of quartets of quantum numbers p and q. The 
spin-less density function (01 calculated from p(r) = ( v F|5(r)| v F), using the second quantized 
form of the operator © 

6(r) = yy p a q 6 {^ p (r')\—L- 6(r - /) 6(0 - 6(<p - <p'M q (r')) 
pq 

= J]ala q S^^^rW^W.vyR^W^W.ip) , (17) 



yields 



P(r) = Ya D » 6 V« K^r^ ^yR^Yi^,^ . (18) 



To illustrate the spherical harmonics content of the density in the Hartree-Fock approxima- 
tion, consider the atomic term \s 2 2p 1 ( 3 P)3d A F for which the {Mi, Ms) — (+3, +3/2) subspace 
reduces to a single Slater determinant 

y a iSM L M s = ®(ls 2 2p 2 ( 3 P)3d 4 F +X+V2 ) = \lsTs2p +1 2p 3d +2 \ . (19) 

When evaluating ( fTSI ). all non-zero Devalues appear on the diagonal (p = q), yielding 

p(r) 4 ^^ = \^ ls (r)\ 2 + \^r)\ 2 + |^ 2p+1 (r)| 2 + |^ (r)| 2 + |^ +2 (r)| 2 . (20) 



This density has a clear non-spherical angular dependence. However, referring to 1 20], 

j j 

W* M m = \Yjm(A <p)f = b "(J' M ) P 2n(C0S ff) = KV> M ) y 2„ o(#, f) (21) 
n=0 n=0 

one recovers the even Legendre polynomial content of the density, although not reaching the 
(2L = 6) limit T6o(#> <p) of the exact density ( TT3l . However this limit will be attained when 
extending the one-electron orbital active set to higher angular momentum values for building a 
correlated wave function. 



Mixed contributions (p + q) may appear in (fTSI ) through ofF-diagonal matrix elements in 
the CSF basis. For example, the interaction of <&(\s 1 2p 1 { 3 P)3d -F+3,+3/2) with the angular 
correlation component <&(ls 2 2p3d( i F)4f 4 F+3,+3/2) , a single electron excitation 2p — > 4/, 
gives rise to Y* Q Y 3[j and Y! +1 Y~ + , contributions. But these contributions are also limited to even 
Legendre polynomials, as appearing in equation (13[ . Indeed, starting from the Clebsch-Gordan 
series [20] 



h+h 1 

Y hmi (&,<p)Y km2 (V,<p)= Yj Z 
l=\h-h\m=-l 



(2h + l)(2l 2 + 1)(2/+ 1) 



4-n 



1/2 


h h I 


1 


(-1)'" 




V ) 


V 



h h 



I 



mi mi —m 
(22) 



and using 



y,_ m (i? ) ^ = (-i) m y,;(i?, V ), 



(23) 



one finds that any contribution of the type Y* Y^ arising from a single electron excitation \liq) 
\hq) preserving the parity, ie. (-1)' 1 = (-1)' 2 , takes the form 



h+h 



f; i? (i?,^ 29 (^,0 = (-i)« J] 

i ac «=\n-h 



(2h + l)(2l 2 + 1)(2Z + 1) 



4n 



1/2 



h h I 




h 



h I 



-q +q 
(24) 



At this stage, we would like to stress that in an MCHF calculation the density never contains - 
what Fertig and Kohn I17I called - "offending" spherical harmonic components, whatever the 
maximum /-value of the orbital active space. 

2.3. The spherical density function 

A spherically symmetric density function can be obtained for an arbitrary CSF $> a Ls m l m s by 
averaging the (2L + l)(2S + 1) magnetic components of the spin-less density function 

1 



P(r)" = 



(2L+ 1)(25 + 1) 



2>« 



LSM L M S 



where p(r) LSMlMs is constructed according to eq. ( fT8l 

p(j) lsm l m s = J^{® aLSMLMs \ala q \® aL s MLMs ) 6„ lspMsq i^(r)i^(r) . 



(25) 



(26) 



Applying equations d25l l and (|26T i for the atomic term ls 2 2p 2 ( 3 P)3d 4 F considered in the 
previous section, we simply get 



p(r)4F = i i 2p2 ^ r) + 2P i (r) + p2 M 



(27) 



which is, in contrast to eq. d20b . obviously spherically symmetric. The sum over (Ml, Ms) 
performed in (l25l > guarantees, for any nZ-subshell, the presence of all necessary components 
{F/ m , | mi = —/,... + /} with the same weight factor, which permits the application of Unsold's 
theorem I21I 

, 21 + 1 

J] \Y lmi (§,<f)\ 2 = — — (28) 



4n 

mi=-l 



and yields the spherical symmetry. This result is valid for any single CSF 

1 

Anr 2 



^ a = iZ«"^w- (29) 

where q„\ is the occupation number of nZ-subshell. Its sphericity explicitly appears by rewriting 

1 as 

p(r) = p(r) \YaoW, iff = ^ |T o(^, 0| 2 , (30) 

p( r ) = ^Yj q "' p2 -' (r) ^ (31) 

«/ 

D(r) = r 2 p(r) = ^ q„iP 2 nl (r) = J] <?„, r 2 ^2 ; (r) . (32) 



with 



and 



The radial distribution function D(r) represents the probability of finding an electron between 
the distances r and r + dr from the nucleus, regardless of directioij^. This radial density function 
reveals the atomic shell structure when plotted as function of r. Its integration over r gives the 
total number of electrons of the system 

D(r) dr = r 2 p(r) dr =Yq nl = N. (33) 
Jo V 

Where above the spherical symmetry of the average density ( f25l ) is demonstrated for a single 
CSF thanks to Unsold's theorem, it can be demonstrated in the general case by combining (1251 ). 
O and the 3-;' sum rule Q 



I<-» 



L-M L 



L k L 



-M L M L J 



(2k+\) ll2 5 kfi (34) 



2 Note that, although denoted as D, this function (evaluated at the r = 
density" used in the context of isotope shifts 1 22]. The latter is indeed p(0) = 

9 



0) is not the so-called "modified electron 
4np(0). 



for each k — 21 contribution (fT4l >. However, the radial density p(r) will be more complicated 
than OTT l. involving mixed contributions of the type P n >i(r)P n i(r) = r 2 R n >i(r)R n i(r), as developed 
below. 

Instead of obtaining a spherically symmetric density function by averaging the magnetic 
components p(r) LLSMiMs through eq. d25l ). one can build a radial density operator associated to 
the function d32l l which is spin- and angular-independent, i.e. independent of the spin (cr) and 
angular {&, <p) variables. Adopting the methodology used by Helgaker et al [13] for defining the 
spin-less density operator, we write a general first quantization spin-free radial operator 

N 

f = J]f(n) (35) 

!=1 

in second quantization as 

pi 

where f pq is the one-electron integral 



f pq = j i/r* (x)/(r)(/r 9 (x)r 2 sin MrdMipdcr . (37) 



Applying this formalism to the radial density operator 

N 

5(r) = 25(r-r,), (38) 

i=i 

and using the spin-orbital factorization (|2|i for both p and q quartets, we obtain the second quan- 
tization form 

$(r) = Y J d Pq (r)ala q , (39) 
pi 

with 

d P q(r) = 6 Vq 6 mipmiq S mgpm3q Rl nln (r)R„ q i q (r)r 2 , (40) 

where the Kronecker delta arises from the orthonormality property of the spherical harmonics 
and spin functions. With real radial one-electron functions, the operator (l39l becomes 

<5(r) = ^ Si>i 6 m > m! 5 m>i <£ n , Vlrim fnimxm, R„'i'(r)R„i(r)r 2 (41) 

'lmim s n ^ m / m s 

R„-i(r)R n i(ry . (42) 
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Its expectation value provides the radial density function D(r) = r 2 p(r) = 4nr 2 p(r) defined by 
ED and (02). 

Building the coupled tensor of ranks (00) from the [2(21 +1)] components of the creation and 
annihilation operators I23I 



\(00) 

/oo 



V2(2/+ 1) 



lm,m s anlm i m s 



I f V 

the operator (l4TT i becomes 



(43) 



(44) 



The expectation value of this operator provides the spherical density function for any atomic 
state. Note that, in contrast to (12611 . the tensorial ranks (00) garantee the diagonal character in 
L, S, Ml and Ms , thanks to Wigner-Eckart theorem 



(aLSM L M s \T^\a'L'S'M' L M s ) = (-1) 



L+S-M,-M s 



L U 
-M L M' L j 



SOS' 
-M s M' s j 



Moreover, the Ml/Ms independence emerges from the special 3 j'-symbol 



j / 



-m, m' } 



{aLS\\T m \\a'L'S'} 
(45) 

(46) 



In other words, where the non-spherical components are washed out by the averaging process 
(1251 1. they simply do not exist and will never appear for the density calculated from (1441 . for any 
(Ml, Ms) magnetic component. 

The radial distribution function D(r) = r 2 p(r) can be calculated from the expectation value 
of the operator d44l >. using the wave function (Q3 or In the most general case (expansion (f3])), 
using the (L5)7-coupled form of the excitation operator, 



fa 1 " a ,f 0)0 - a ,Y° 0) 



(47) 



one obtains 



WamWWam) = (-1) 7 



7 7 
-MOM 



(48) 
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with 

F*j» = -£ V2(2TTTj X (a>„/f° /, («'/, nl) , (49) 

/= 1 n,n' 

and 

I p (n'l,nl)(r) = R^R^r)/ 2 . (50) 

The diagonal reduced matrix element (RME) evanuated with the Breit-Pauli eigenvector (O has 
the following form 

(VojWF^Woj) = J] afo <<D(a,L i 5 ; 7)||^ 00)0 ||<D(a i L j 5 i y)) (51) 

U 

where the RME in the (LS)J coupled basis reduces to 



maiUSiJMW^majLjSjJM))^ ^^—^ll— ^ {<b( ai LiS dW^majLjS j))6 Ll , Li 5 Sl ,s s 

(52) 

and 

^ = -2 V^Tlj J] (a^a^C 7 P (n'Z, nZ) . (53) 

l=\ n,ri 

From the analogy of the operator (l53l and the non-relativistic one-body Hamiltonian op- 



erator (see eq. (A5) of J^IX one observes that the angular coefficients of the radial functions 
I p (n'l, nt) (r) are identical to those of the one-electron Hamiltonian radial integrals /„</,„/ , as an- 
ticipated from McWeeny analysis [3]. These angular coefficients can be derived by working out 
the matrix elements of a one-particle scalar operator F p between configuration state functions 
with u open shells, as explicitly derived by Gaigalas et al 12511 who expressed them as a sum over 
one-electron contributions 

{<S>(aLS) ||^ 00) | ®(a'LS)) = J] <<D(aLS) l^fo/,-, njlAl <S>(a'LS)) (54) 



where 



(<t>(aLS)\\F p (nl i ,n j lj)\\®(a'LS)} 

= (-1) A+I V2(2/,- + 1) R (A h Aj, A hra , A ket ) 6 Uj I p (mh, n jlj) 

+(1 - 8{n.i,nj)) Uiilf aiQjLjSi ^a^ 5 ^ n ^i' a iQi L iS,) 

x {n/]' ajQjLjSj \\a^\\ n j]' ajQjLjS ,)} . (55) 
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In this last expression, A = I or s, (<f>(aLS)\ and \(S>(a'LS)) are respectively bra and ket functions 
with u open subshells, 

K bm = (LiSuLjS j,Li,Si,,LyS j,) bm and A ke ' = (i^S h LjS j,Li>S ?,L r S ff" denote the respec- 
tive sets of active subshell angular momenta. The operators o„ are second quantization op- 
erators in quasispin space of rank q — 1/2. The operator c^^ m m — dltfm, creates electrons 
with angular momentum quantum numbers l,mi,s,m s and its conjugate afi^mim ~ ^mm s = 
(—l) l+s ~ mi ~ ms a ( Lm,m annihilates electrons with the same quantum numbers I, mi, s,m s in a given 
subshell. The coefficient R (Aj, Aj, A bra , A to ) is the recoupling matrices in /- and s- spaces and A 
is a phase factor. 

3. Density matrix and natural orbitals 

Using 148) . ( BIT ). ( 152b and ( 1331 ). the radial distribution function gets the following form 



D(r) = r 2 p(r) = £ a*DyW«j = J] fl <* Z Z v ™'' / p ( "' / ' "° ^ ' (56) 

ij ij Y I n'n 

which can be rewritten in a compact form 

D ( r )=ZZ jo »'» / >' u/) ' (57) 

/ n'n 

with 

Pn>n = Yj a *i V nn>l a j ■ {5 %) 
U 

The Si, j. Kronecker appearing in 05] ) assures the block-structure of the density matrix p whose 
elements are defined by ( f58l > for the /-angular symmetry. 

The natural orbitals (NO) are defined as the one-electron functions that diagonalize the den- 
sity matrix p 

C 'p C = p . (59) 
Within a specific angular /-symmetry, the eigenvalue problem for the relevant /-block 

p'C 1 = C'p 1 (60) 

defines the natural radial orbitals through the following transformation 

n 

The eigenvalues [A 1 , = pL] are interpreted as the occupation numbers of the NOs {Rki(r)}. 
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4. Algorithm description 

To calculate the radial density function and the natural orbitals from an arbitrary A^-electron 
wavefunction ^ajM, we wrote a FORTRAN implementation of equation (|48| |. as an extension of 
the atsp2K package. The essential part in the calculation of the density function, is the evaluation 
of the reduced matrix element (f5Tb . In pseudo-code, the reduced matrix element dSTb is written 

as 



< VaJ II ^ 00) ° II VaJ > = J] Z a>aj Z Z 7 >' v) UNITELEMENT^, v) 

i j 11 v 

SPIN_ANGULAR_DENSITY(CSF,, //; CSF y , v; 00) . (62) 

where 

UNITELEMENTOu, v) = -[^, s M ]*<S(Z M , ? v ) . (63) 

The routine SPIN_ANGULAR_DENSITY, is inspired by the routine NONHIPER of the hf s hy- 
perfine structures program of ATSP2K. It organizes the calculation of the spin-angular part of 



ONEPARTI- 



by calling the subroutine ONEPARTICLE 1 or ONEPARTICLE2 from 
CLE1 performs the calculation of the spin-angular part when the one-electron operator acts on 
one open shell and ONEPARTICLE2 performs the calculation when the operator acts on two 
open shells. Both calculate the spin-angular part using the expresion (l55l l in which /p(n,/,-, njlj) = 
1 . The products of the weight factors with the corresponding spin-angular part are stored and 
accumulated in the two dimensional array FACTORMATRIX(//, v) where the rows and columns 
are defined by the {nl) subshell quantum numbers of the bra and ket, respectively. FACTORMA- 
TRIX is the precursor of the density matrix (|58V The products of the array elements with their 
corresponding radial part I p (n'l, nl) are accumulated to build the radial distribution function ( f57l i. 
The reader is referred to the flowchart in figure Q] for a schematic overview of the calculation of 
the density function. 

The NOs are obtained by diagonalizing this matrix and using the eigenvectors to construct the 
orbitals. The diagonalization of the density matrix d59l l is performed using the DSYEV subroutine 



from the Lapack [26] library. This routine computes all eigenvalues and eigenvectors for a given 
real symmetric matrix. The NOs are ordered and labelled according to their occupation numbers 
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Most of the subroutines needed for density exist in hf s of ATSP2K, besides the routines 
from the ATSP2K libraries. The new modules are density . f , spin_angular .density . f and 
unitelement . f . The code readwfn.f that reads in the wave functions differs from the one 
encountered in hf s by the COMMON/ ADATA2/ AT , TT , ELNAME (NWD) needed to store the ATOM, 
TERM and ELNAME variables. 

As an illustration, an interactive session is described in appendixlAl for an — 3 CAS-MCHF 
expansion of the beryllium ground state (63 CSFs). Upon execution of density, the user is 
asked to specify the name of the data files, which were obtained from an ATSP2K run. density 
then reads the CSF weights ({c,} and {a,} for the non-relativisic and Breit-Pauli expansions, re- 
spectively), the configuration state functions quantum numbers and the radial functions from the 
files. The conventions of the data and the file types, summarized in table [U were adopted from 
ATSP2K. However, for a relativistic calculation, the . j should be renamed file . 1 and edited to 
extract the selected relativistic /-eigenvector of interest. 

In an interactive session, density asks the user a few questions concerning the output and 
wether the NOs should be evaluated. In table|2]we list and comment the questions. Most of the 
output, however, is written to disk. The output files produced by the program are summarized 
in table [3] The name . d file, which is always generated, contains the radial distribution D(r) and 
density p(r) functions. The density program by default generates some output to the standard 
out: the "modified electron density" [22] at the nucleus (p(0) = 47rp(0)), the occupation numbers 
of the natural orbitals, with their composition in terms of the original orbitals, and as a final 
check, the integral of the density function that should give the total number of electrons accord- 
ing to d33l . The P„i(p) = r~ l t 2 P„i(r) functions appearing in the files name. pit and name.n 
are defined in the logaritmic variable p = log e (Zr) 111 ill for the original and natural orbitals re- 
spectively. If the user asks for more details ( ' yes ' to the question PRINT ALL DATA (y/*) ) , 
density prints out the contributions to the reduced matrix element d62l . providing for each pair 
(z, J) of CSFs, the labels (p, v) of the orbitals involved, the corresponding spin-angular coeffi- 
cient, together with the relevant weights product (a,fl)). Using this option, the user also gets the 
contributions to the modified density at the nucleus, the norm of the input orbitals, the matrix 
elements of the density matrix, and the natural orbitals (before they are sorted according to their 
occupation number), with their complete eigenvector composition. 

To install the program (FORTRAN 90 compilation and linking with the ATSP2K libraries), 
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the provided Install script should be edited to set the appropriate path and environment variables. 



5. Applications and examples 

To illustrate the data in the output files, we plotted in figure [2] the radial density distribution 
D(r) = r 2 p(r) from the . d output file calculated for a CAS-MCHF wave function of the beryllium 
ground state (Be ls 2 2s 2 l S), using a n = 9 orbital active set. In the same figure, the Hartree-Fock 
radial density is compared with the one obtained with two correlation models: i) the n — 2 CAS- 
MCHF expansion, largely dominated by the near-degeneracy mixing associated to the Layzer 
complex ls 2 {2s 2 + 2p 2 } and ii) the n — 9 CAS-MCHF. From the plotted results we notice that 
the density of the n — 2 calculation already contains the major correlation effects, compared to 
the n-9 calculation. Indeed, the density does not seem to change a lot by going from the n - 2 
to the n-9 orbital basis, the valence double excitation \ s 2 2p 2 contributing for 9.7% of the wave 
function. From the energy point of view however, this observation is somewhat surprising (see 
table 0J: the correlation energy associated to the n = 2 CAS-MCHF solution "only" represents 
47% of the n-9 correlation energy. 

In a separated pair-MCHF approach, the reduced forms of the CSF expansions are often used 
to get a compact multiconfiguration representation of the state and to avoid possible variational 
redundancies between orbital rotations and mixing coefficients transformations. For some spe- 
cific cases, the so-produced MCHF one-electron functions are nothing else than the natural or- 
bitals . For expansions closed under orbital rotations, one can test our density computational 
tool by: 1) perfoming an (unreduced) MCHF calculation, 2) obtain the natural orbitals from the 
diagonalization of the density matrix and 3) making a CI calculation in the resulting NO basis. 
Both calculations should yield the same total energy for two rather different representations of 
the same total wave function. Amongst the two, the NO-CSF expansion is naturally condensed. 
This is illustrated in table|5]for a n = 5 SD-MCHF valence correlation calculation on the ground 
state of Be (E = -14.619 083 a.u., using a Hartree-Fock frozen core). The eigenvectors calcu- 
lated in both MCHF and NO one-electron bases are reported and compared to each other. Note 
that, in this specific case (a pair of S e symmetry), the transformation that diagonalizes the den- 
sity matrix eliminates the off-diagonal (n + n') contributions \s 2 nln'l [27]. The reduction in the 
number of CSFs (30 — > 15) through the use of NOs is quite impressive. For a n — 6 SD-MCHF 
valence correlation calculation the CI-NO approach yields a CSF expansion with 29 terms less 
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and for a CAS-MCHF n = 9 wave function (271 733 CSFs), the NO basis leads to a reduction of 
15 695 CSFs. 

As a third example, we illustrate the influence of relativistic effects - in the Breit-Pauli ap- 
proximation - on the density function of the Be-like +4 atom, by comparing the densities of 
the fine-structure states \s 2 2s2p 3 Pq, 3 Pj and 3 P° 2 . From the plots in figure[3]and the data given 
in table [6] we observe that the largest energy difference corresponds to the largest difference in 
density function. More bound is the level, higher is the electron density in the inner region. 

When studying the electron affinities, it is often interesting to investigate the differential 



correlation effects between the negative ion and the neutral system 12811 . Figure [4] displays the 
density functions D(r) of both the [Ne]3s 2 3/? 4 3 P ground state of neutral Sulphur (S) and the 
[Ne]3s 2 3/? 5 2 P° ground state of the negative ion S~, evaluated with elaborate correlation models 
B29I1 . together with their difference AD(r). The latter reveals where the "extra" electron lies and 
its integration gives one, as it should. 
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A. An interactive session 

Scat n3.c 



ls( 2) 2s( 2) 
ISO ISO IS 

ls( 2) 2s( 1) 3s( 1) 
ISO 2S1 2S1 2S IS 

ls( 2) 2p( 2) 
ISO ISO IS 

ls( 2) 2p( 1) 3p( 1) 
ISO 2P1 2P1 2P IS 

1st 2) 3s( 2) 
ISO ISO IS 

ls( 2) 3p( 2) 
ISO ISO IS 

ls( 2) 3d( 2) 
ISO ISO IS 

ls( 1) 2s( 2) 3s( 1) 

3p( 4) 
ISO 

3p( 2) 3d( 2) 
ISO ISO IS 

3p( 2) 3d( 2) 
1D2 1D2 IS 

3p( 2) 3d( 2) 
3P2 3P2 IS 

3d( 4) 
ISO 

3d( 4) 
1S4 



$cat n3.1 

Be Z = 4 . NEL = NCFG = 63 



2*J = NUMBER = 1 
Ssms = 0.484179758 
1 -14.654414586 ls(2) .2s(2)_lS 
0.95181933 0.00029779 0.30027819 0.00037936-0.00118903-0.00023749-0.01763502 
0.00019391-0.04365498 0.00316223-0.00663489 0.00377678-0.00043175 0.00173694 
-0.00032836-0.00086628 0.00002903 0.00159570-0.00108079-0.00181964 0.00002085 
-0.00001498-0.00000941 0.00415462 0.00703251-0.02349037 0.02848932-0.00012045 
-0.00265663 0.00007304-0.00016224-0.00019179 0.00009463 0.00016368 0.00001426 
0.00011178-0.00003997 0.00061770 0.00194285-0.00725367-0.00004543 0.00897260 
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. 00000543-0 . 00001468-0 . 00005233-0 . 00000271-0 . 00000675-0 . 00000998-0 . 00002343 
. 00003445 . 00000814-0 . 00013590-0 . 00000014-0 . 00000679-0 . 00002475 . 00042840 
0.00000708-0.00000844-0.00052396 0.00000473 0.00000004 0.00000119 0.00000002 

$density 

Density calculation, Summer 2009 

Give <name> of the <name>.c, <name> . 1 <name>.w files: 
n3 

Files: n3 



PRINT THE ORBITALS (*/n) 

Printout orbitals 

PRINT THE MATRIX O/n) 

Printout the matrix 

CALCULATE NATURAL ORBITALS O/n) 

Calculate natural orbitals 

PRINT ALL DATA (y/*) 

Do not print all informations 

ANALYSING THE CALCULATION 



ACCURACY IS SET TO 1 . 0000000000000007E-016 



STATE (WITH 63 CONFIGURATIONS) : 



THERE ARE 6 ORBITALS AS FOLLOWS: 

Is 2s 2p 3s 3p 3d 
THERE ARE CLOSED SUBSHELLS COMMON TO ALL CONFIGURATIONS AS FOLLOWS: 



NORM OF WEIGHTS = 1.000000004562740 
ATOM Be TERM lSe 

ALL WAVEFUNCTIONS EXIST. 



START OF THE DENSITY CALCULATION 



MODIFIED ELECTRON DENSITY AT THE NUCLEUS: 
= 444.31734212383130000 

EIGENVECTOR: 

1 = Eigenvalue 6 : . 19968595313710157E+01 
Is ' = 

-0.99450714610441153E+00 Is AZ= . 14887071657598840E+02 

0.10466865508139691E+00 2s AZ= . 10194455194727872E+01 

-0.94819357413400507E-04 3s AZ= . 23772474518338814E+02 

2 = Eigenvalue 5 : . 18147702141513149E+01 
2s ' = 

0.99450712354440736E+00 2s AZ= . 10194455194727872E+01 
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0.10466862950576727E+00 Is AZ= . 14887071657598840E+02 

0.24334504925053317E-03 3s AZ= . 23772474518338814E+02 

3 = Eigenvalue 2 : . 1250344482753356BE-02 
3s ' = 

0.99999996589623770E+00 3s AZ= . 23772474518338814E+02 

-0.11976912542329965E-03 Is AZ= . 14887071657598840E+02 

-0. 23208377825771 107E-03 2s AZ= . 10194455194727872E+01 

4 = Eigenvalue 4 : . 18458626963217419E+00 
2p ' = 

-0.99999853877956664E+00 2p AZ= . 15057313981228271E+01 

-0.17095141798808027E-02 3p AZ= 0.51881186928943286E+02 

5 = Eigenvalue 3 : . 18993473382529018E-02 
3p ' = 

-0.99999853877956664E+00 3p AZ= 0.51881186928943286E+02 

0.17095141798808027E-02 2p AZ= . 15057313981228271E+01 

6 = Eigenvalue 1 : . 63431 127B44789989E-03 
3d ' = 

0.10000000000000000E+01 3d AZ= . 31738718272621771E+00 

SUM OF EIGENVALUES 4.000000018250959 

INTEGRAL OF THE DENSITY FUNCTION: 
N = 4.00000001825096200 

DENSITY FUNCTION IS IN FILE n3.d 
END. 
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extension data in the file 



. c configuration state function (CSF) expansion 

. w radial wave functions (numerical values in binary form) 

. 1 expansion coefficients from a non-relativistic (LS) calculation 

. j expansion coefficients from a Breit-Pauli (LS J) calculation 



Table 1 : File convention 



Answer Implication 



Question 

PRINT THE ORBITALS (*/n) y 

PRINT THE MATRIX (*/n) y 

CALCULATE NATURAL ORBITALS (*/n) y 

PRINT ALL DATA (y/*) y 



The input radial functions will be written to . pit . 

The density matrix will be written to . matrix . 

Calculate the NOs and write them on . n (formatted) 

and . nw (unformatted) files. 

Detailed output written to std out: 

MODIFIED DENSITY AT THE NUCLEUS 

NORM OF THE ORBITALS 

DENSITY MATRIX 

EIGENVALUES AND EIGENVECTORS 



Table 2: Questions density asks the user. "*" indicates the default answer. 



extension data in the file 

Tplt n, Rniin), P n i(ri) = rtRnin), P nl {pi) = r- lll P nl (n) 

.d r i ,p(r i ),D(r i ) = rfp(r i ) 

. n n, R n i(n), Pni(n) = nR^n), Pniipd for Natural Orbitals 

. nw analogue of . w for the Natural Orbitals (contains P„/(p,)) 



Table 3: Output files created by density 



model energy (a.u.) correlation energy (a.u.) 

HF -14.573 023 

n = 2-CAS -14.616 856 E n = 2 - E HF = 0.043 832 

«=9-CAS -14.667 013 E n = 9 - E HF = 0.093 986 



Table 4: Total energy for the ground state of Be with different correlation models. 



21 





CSF 




ML.Hr basis 


Natural orbital basis 


ls( 


2) 2s( 2) 




0.95282855 


0.95370264 


ISO 


ISO IS 








ls( 


2) 2s( 1) 3s( 


1) 


0.03858929 


0.00000000 


ISO 


2S1 2S1 2S IS 








ls( 


2) 2s( 1) 4s( 


1) 


-0.01524193 


0.00000000 


ISO 


2S1 2S1 2S IS 








ls( 


2) 2s( 1) 5s( 


1) 


0.00133387 


-0.00000001 


ISO 


2S1 2S1 2S IS 








ls( 


2) 2p( 2) 




0.00133387 


0.29736974 


ISO 


ISO IS 








ls( 


2) 2p( 1) 3p( 


1) 


-0.00032489 


0.00000000 


ISO 


2P1 2P1 2P IS 








ls( 


2) 2p( 1) 4p( 


1) 


-0.00019862 


0.00000000 


ISO 


2P1 2P1 2P IS 








ls( 


2) 2p( 1) 5p( 


1) 


0.00089172 


-0.00000001 


ISO 


2P1 2P1 2P IS 








ls( 


2) 3s( 2) 




-0.03930620 


-0.04031077 


ISO 


ISO IS 








ls( 


2) 3s( 1) 4s( 


1) 


-0.00463218 


0.00000000 


ISO 


2S1 2S1 2S IS 








ls( 


2) 3s( 1) 5s( 


1) 


0.00091733 


0.00000000 


ISO 


2S1 2S1 2S IS 








ls( 


2) 3p( 2) 




0.29736945 


0.00532117 


ISO 


ISO IS 








ls( 


2) 3p( 1) 4p( 


1) 


-0.00003969 


0.00000000 


ISO 


2P1 2P1 2P IS 








ls( 


2) 3p( 1) 5p( 


1) 


0.00024549 


0.00000000 


ISO 


2P1 2P1 2P IS 








ls( 


2) 3d( 2) 




-0.01669194 


-0.01669247 


ISO 


ISO IS 








ls( 


2) 3d( 1) 4d( 


1) 


0.00005217 


0.00000000 


ISO 


2D1 2D1 2D IS 








ls( 


2) 3d( 1) 5d( 


1) 


-0.0001 1379 


0.00000000 


ISO 


2D1 2D1 2D IS 








ls( 


2) 4s( 2) 




-0.00422002 


-0.00432946 


ISO 


ISO IS 








ls( 


2) 4s( 1) 5s( 


1) 


-0.00107435 


0.00000000 


ISO 


2S1 2S1 2S IS 








ls( 


2) 4p( 2) 




0.00182955 


0.00184355 


ISO 


ISO IS 








ls( 


2) 4p( 1) 5p( 


1) 


0.00032865 


0.00000000 


ISO 


2P1 2P1 2P IS 








ls( 


2) 4d( 2) 




-0.00361419 


-0.00363174 


ISO 


ISO IS 








ls( 


2) 4d( 1) 5d( 


1) 


0.00030308 


0.00000000 


ISO 


2D1 2D1 2D IS 








ls( 


2) 4f( 2) 




0.00618640 


0.00621375 












ls( 


2) 4f( 1) 5f( 


1) 


-0.00048678 


0.00000000 


ISO 


2F1 2F1 2F IS 








ls( 


2) 5s( 2) 




-0.00160546 


-0.00136552 


ISO 


ISO IS 








ls( 


2) 5p( 2) 




-0.00141498 


-0.00149216 


ISO 


ISO IS 








ls( 


2) 5d( 2) 




-0.00103723 


-0.00101914 


ISO 


ISO IS 








ls( 


2) 5f( 2) 




0.00188266 


0.00185530 


ISO 


ISO IS 








ls( 


2) 5g( 2) 




-0.00284386 


-0.00284386 


ISO 
* 


ISO IS 









Table 5: Comparison of Be n = 5-valence eigenvectors in the MCHF and NO bases. 
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model energy (a.u.) energy difference (a.u.) 

Is z 2s2p J P° -68.032 086 

3 P° -68.031 473 AE W = 0.000 613 

3 P° 7 -68.030 102 AE 2 i = 0.001 370 



Table 6: Fine structure total energies of +4 ls 2 2s2p 3 P° 



Read and parse the configuration state function name . c 
Read and parse the mixing coefficients name . 1 
Read and parse the wavefunction name . w 






For each CSF in the Bra ( x ¥ a jM\ 






For each CSF in the Bra |Y ff/M ) 






Calculate the weight product a,a ; 
and the spin-angular part. 
Store the product in the FACTORMATRIX, 
entry («,/,, rtjlj). 






For each entry in the FACTORMATRIX 




Add to the density the product of 
the factor matrix entry times the 
radial functions I p {ji, v) 


Write density function to name . d 



Figure 1 : Flowchart of the Density program 
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r/a 

Figure 3: Comparison of the ls 2 2s2p 3 PJJ, 3 Pj and 3 Pj radial density functions of +4 . Density differences have been 
scaled by a factor 10 000. 
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Figure 4: Ground state S and S density functions | 29]. Density differences have been scaled by a factor 30. 
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